%% fig7a_main.m
% 在 (m21, m12) 参数平面上计算系统的最大 Lyapunov 指数，并绘制彩色热力图
% 固定 FN 神经元参数 ε = 4.87
% 使用封装函数 LEs 计算最大 Lyapunov 指数

clear; clc; close all;

%% 1. 设置固定的系统参数
% HR 神经元1 参数
params.a1 = 1.5;    % 激活参数
params.b1 = 3;      % 非线性抑制系数
params.c1 = 1;      % 恢复变量
params.d1 = 5;      % 抑制作用
params.i1 = 0.45;   % 外部电流

% FN 神经元 参数（FitzHugh-Nagumo 模型）
params.a2 = 0.77;
params.b2 = 0.3;
params.c2 = 0.8;
params.i2 = 0;

% HR 神经元2 参数
params.a3 = 0.9;
params.b3 = 3;
params.c3 = 1;
params.d3 = 5;
params.i3 = 0.5;

% 固定耦合系数（其他耦合参数）
params.m23 = 0.785; % FN -> HR2
params.m32 = 1.01;  % HR2 -> FN

% 注意：接下来扫描的参数为 m21 和 m12

%% 2. 定义参数扫描范围
m21_vals = linspace(2.5, 3.5, 100);   % m21 取值范围
m12_vals = linspace(4, 5, 100);         % m12 取值范围

%% 3. 设定积分与 Lyapunov 计算参数
T_trans = 1000;   % 去除瞬态时间
T_le = 100;       % 用于计算 Lyapunov 指数的时间
dt_save = 10;     % 每隔 dt_save 时间正交化一次

RelTol = 1e-6;    % 相对容忍度
AbsTol = 1e-8;    % 绝对容忍度

% 固定的 FN 神经元参数 ε（论文中给定）
eps_val = 4.87;

% 初始条件（与论文一致）
X0 = [-2, 0, 0, 0, 0, 0.1, 0.3];

%% 4. 预分配存储最大 Lyapunov 指数的矩阵
nm21 = length(m21_vals);
nm12 = length(m12_vals);
MLE_map = zeros(nm12, nm21);  % 行对应 m12, 列对应 m21

%% 5. 并行计算 (m21, m12) 网格上的最大 Lyapunov 指数
if isempty(gcp('nocreate'))
    parpool;
end

parfor i = 1:nm21
    m21_cur = m21_vals(i);
    rowData = zeros(1, nm12);
    for j = 1:nm12
        m12_cur = m12_vals(j);
        % 在 parfor 中避免直接修改共享变量，创建局部副本
        localParams = params;
        localParams.m21 = m21_cur;
        localParams.m12 = m12_cur;
        
        % 调用封装函数 LEs 计算前 numLE 个 Lyapunov 指数
        LE_val = LEs(eps_val, localParams, X0, T_trans, T_le, dt_save, RelTol, AbsTol, 2);
        rowData(j) = LE_val(1);  % 只记录最大 Lyapunov 指数
    end
    MLE_map(:, i) = rowData;
end

%% 6. 绘图
figure;
imagesc(m21_vals, m12_vals, MLE_map);
axis xy;           % 确保 y 轴向上
colorbar;          % 显示色标
colormap jet;      % 使用 jet 色图
xlabel('m_{21}', 'FontSize', 12);
ylabel('m_{12}', 'FontSize', 12);
title('Fig.7(a): Maximum Lyapunov Exponent in (m_{21}, m_{12}) Plane', 'FontSize', 14);
set(gca, 'FontSize', 12);
